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I. INTRODUCTION 


Thermal stresses have become increasingly important in 
engineering practice during recent years. In power genera- 
tion higher cycle temperatures and use of nuclear fission 
are largely responsible for this trend. This thesis describes 
a computer program for finding temperatures and stresses in 
bodies having axisymmetric geometry and loading. The govern- 
ing equations are those of isotropic, uncoupled, quasi-static, 
linear thermoelasticity. They are discretized by using the 
finite element method. A FORTRAN IV computer program using 
double-precision arithmetic has been written to solve 


problems of the following kinds. 


A. TEMPERATURE PROBLEMS 

The transient temperature vector, evaluated at the nodal 
points, may be obtained for an axisymmetric body with the 
combination of insulated, convection, or constant temperature 
boundary conditions. For the convection thermal boundary 
condition, however, we may have fluid flowing with entry 
temperature prescribed as a linear function of time (RAMP). 
The program can handle up to 15 different ramps, each having 
a different flow velocity, and with discontinuities between 


successive ramps. 


B, STRESS PROBLEMS 
The program will generate load vectors for pressure load- 


ing, centrifugal loading, and axial force, Provision is 


Sern renner 


made for direct input of one additional load vector. Stresses 


may be found for any combination of these loadings. 


C. THERMAL STRESS PROBLEMS 

* Thermal stresses may be found for as many as 20 different 
temperature vectors which may be output of the temperature 
problem or direct input. In short, in this pare any combina- 


tion of the temperature and stress problems may be used. 


II. . FINITE ELEMENT FORMULATION OF HEAT CONDUCTION 
IN AXISYMMETRIC BODIES 


A. METHOD OF FORMULATION 
For bodies o£ revolution under axisymmetric loading the 
mathematical problems presented are two-dimensional. The 


governing equation for non-steady heat conduction is 
V-kVT = oct, (1) 


where k is the thermal conductivity, p the density, c the 
specific heat, T the temperature and V the gradient operator. 
The superior dot denotes a time derivative. 


Applying Galerkin's principle [1] gives 


s N,V-KVT dv = Spe Ny T dV, (2) 
where the integral is over the volume V of the conducting 
body and N; is a "shape function" used in representing the 
temperature distribution. If S is the surface of the body 
and n the outward normal to surface, then using Gauss' 


theorem we can write 


f Ve (N. kT) aV= Ss N. Ket ds (3) 
V * S u 
Since 
£¥°(N; VT) dv = sf (WN,)* (KVT) dv 
V Vv 
+ S NW (kVT) av, (4) 
V 1 


we can combine Eqs. 2, 3 and 4 to get 


; spas tee eet oT 
f pcN, T av + f (WN) (KPT) av = J Ni ks— ds. (5) 


Each node of the solid region has a separate discretized 
linear equation calculated from Eq. 5 using the appropriate 
shape function N,. Thus each of the volume integrals on the 
left hand side of Eq. 5 yields a square coefficient matrix 
in the assembled set of equations. Calculation of these 
matrices is a standard process. Details are given by 
Zienkiewicz [1]. 


The discretized set of equations takes the form 


Ile 
Je 


+YT=v, (6) 


where there is a term by term correspondence with Eq. S. 
The real symmetric matrices C and Y represent, respectively, 
the thermal capacitance and thermal admittance. The elements 
of the vector T are nodal temperatures. The vector v, dis- 
cussed in the following section, depends upon the thermal 
boundary conditions. 

In the present development piecewise constant material 
properties have been assumed, i.e., k,p and c are constant 
within each element, but may vary from element to element. 


Also two-dimensional isoparametric elements are used. Appli- 


cable equations are summarized in Appendix A. 


ae this text a double underline denotes a rectangular 
matrix and single underline denotes a column vector. 


B. THERMAL BOUNDARY CONDITIONS 


Thermal boundary conditions affect only those scalar 
equations of Eq. 6 which correspond to boundary nodes. 
Accordingly, the vector v is sparse. Also, in a single 
problem it is common to have different thermal boundary con- 
ditions on individual portions of the boundary. In what 
follows the subvectors of v (distinguished by individual 
superscripts) which correspond to separate boundary condi- 
tions are treated individually. 

1. Insulated 

It is clear that for insulated boundary conditions 
the subvector vt) of the right hand side of Eq. 6 corres- 
ponding to this boundary condition is zero, since ae = 0. 

2. Convection 

The heat transfer mechanism occurs in the interface 
of the solid and fluid. If the fluid temperature is 6 and 
the heat transfer coefficient is h, then equating heat con- 


ducted away from the surface to the efflux from the solid 


{2} gives 
-k (Z) = h(T-8), (7) 


where the subscript S means that the derivative is evaluated 
at the surface. 


For constant h: 


SN, k5- ds =h f N,(6-T) ds. (8) 
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In what follows the fluid temperature ® is taken to be a 
specified function of position and time. For purposes of 
discretization, the fluid temperature is specified at a 
discrete number of fluid "nodes." If 65 represents the 

fluid temperature at fluid node j, then the fluid temperature 


along the boundary may be represented by 


6-2 N,6,, (9) 


where the N; are one-dimensional forms of the shape functions 
used for the solid. Substitution in Eq. 8 gives 
sn. k 22 as = 5 y*.(0.-T.) (10) 
s i on j = hy ial: INS aed 
where the summation extends over the surface nodes and the 
% 
coefficients Yij are given by 
& 
Yi 


=h NiN; ds. (11) 


j 
Assembling the contributions from Eq. 10, the subvector 


y (2) for the convection boundary condition may be written 


v2) = y%9, (12) 


The contributions -f Vi57j from Eq. 10 are included by 


augmenting the matrix Y (see Eq. 13 below). 
3. Constant Temperature 
If 6 represents the constant temperature desired at 
the wetted surface, then we can use the convection boundary 


condition and replace h by a big number (say 107°} Since h 


is very large, then for thermal equilibrium the temperature 


Re pe te: 
Be Sele rterreremmmmnemmnnmem neem ee 
Bees 


T at the surface will be forced to equai @. So the subvector 
yO) for the portion corresponding to the constant temperature 
. boundary condition can be obtained from Eq. 12. 
Upon the application of these boundary conditions in 
a single problem, the right-hand side vector v will be com- 
bined from the corresponding subvectors and the finite ele- 


ment discretized equation becomes 
CVey ty, (13) 
+ 
where Y =Y+#Y. 
C. EXACT TIME SOLUTION WITH SPATIAL DISCRETIZATION 
We consider only the solution of Eq. 13 for v = constant 


with T= a at time t= 0, Let T. be a particular solution 


(steady state) with qT. = 0 so that 
t= @) 'y. (14) 


For the homogeneous equation 


cise y Ts (15) 


the assumption T = w exp(-At), where w is a vector and A is 


a scalar, yields the form 
yeweagw. (16) 


It is apparent that Eq. 16 defines an eigenvalue problem. 
Let A be the spectral matrix and W the modal matrix with 


normalization according to 


CWFL, (17) 
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where I is the identity matrix of the same order as C. Now 


let 
T = W exp{-At) b (18) 


where b is a constant vector. Substituting this in Eq. 15 


gives 


Y Wexp(-At) b = CWA exp(*At) bd. (19) 


YNrA; (19") 


and this is guaranteed to be satisfied since A and W are 
spectral and modal matrices for the eigenvalue problem of 


Eq. 16 with W 


normalized according to Eq. 17. 

Returning to the original problem (Eq. 13), the com- 
plete solution may be written as 

T = W (8 + exp(-At) b) (20) 
where 
=WB, (21) 
and 8 is a constant vector. Now 8 may be found (using 


Eq. 14) to be 


peaty 


I 


Ve (22) 
Substituting this result into Eq. 20 and using the initial 
condition gives the result 

Ws | (23) 


The general solution of Eq. 13 may thus be written 


as 


EEN YB wll 


T = WQ- exp(-Avy) At wy + W exp(-Aryw ta. (24) 


For purposes of the present program non-zero compo- 
nents of vector v are to be specified as piecewise linear 
functions of time. During each segment of time history of v 
an analytical solution of Eq. 13 is possible in the form of 
a particular solution plus a complementary solution such as 
Eq. 20. At each node the corresponding time variation of 
temperature will consist of a linear part contributed by the 
particular solution and a sum of n terms representing the 
complementary part. Each of these n terms decays exponen- 
tially with a separate time constant. In principle it is a 
Straightforward process to find each particular solution and 
accompanying complementary solution. 

Contemplated problems may typically have from 10 to 
40 segments required to represent the piecewise linear 

: variation of v. The number of body nodes n will be of the 
order of 100 or more. In view of the number of particular 
solutions required, each accompanied by an individual comple- 
mentary solution of the form given by Eq. 20, it was concluded 
that a numerical solution of Eq. 12 would be considerably 
more economical than an analytic one such as that given by 


Eq. 24. 


D. TIME INTEGRATION 

in this section the relative merits of the Runge-Kutta 
and trapezoidal methods of time integration are discussed. 
Since either of these methods will give an exact result if 


the solution is a linear function of time, investigation 


a FER iy errors meneccay ay srniar-annecemnmaneenn ane arene neers or NAO TU AY OTSA PIE RAE SE TENE IER ATEN PIS SET 


is confined to performance on a single scalar equation 
y + dy = G (25) 


whose solution y = Ve exp(-At) is of the same form as the 


components of the complementary solution (Eq. 20 ). 


1. Runge-Kutta Method 
A method introduced by Runge and subsequentiy elabo- 
rated by Heun and Kutta [3] is widely used for the numerical 
solution of first order ordinary differential equations. 
This algorithm prescribes a sequence of calculations for 
+1 2t time Tia 
of Y; and values of y at intermediate and end points of the 


determining the ordinate Ve = T,+AtT in terms 


interval At. The fourth-order form, which requires four 


evaluations of y, gives for Eq. 25 the result 


y; 2 3 4 
fhe nar + Gt”. Ged”, Ged” (26) 


The right-hand side of Eq. 26 represents the first five terms 
of the Taylor expansion of the exact solution 

(Yaa /%3 = exp(-AAt)) so we may conclude that the relative 
error in each time step is less than moduius of the next 
term: (AAt)°/5!. 

In addition to providing the apparent prospect for 
high precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
integration process without requiring additional computa- 

- tionally expensive matrix decompositions. The attractiveness 


of these two features dictated a thorough exploration of the 


potential of this method for the present application. The 


disqualifying defect which emerged after studying a number 
of examples is readily appreciated from examination of 
Table I. For values of AAt less than 0.5 it is apparent 
that Runge-Kutta scheme affords acceptable engineering 
accuracy. However, when the method is appiied to solution 
of Eq. 13 we must deal with a number of A‘s equal to the 
number of nodes (see Eq. 20 ). This number may be greater 
than 100 and the ratio of the largest \ to the smallest may 
easily exceed 1000. Although the solution is dominated by 
the Sone hens of the eigenvectors corresponding to the 
smaller As, it is clear that the solution will be unstable 
if the largest \At exceeds about 2.7. Because an unaccept- 
ably small At is required in typical problems, the Runge- 
Kutta method was rejected. 

2. Trapezoidal Method 


The trapezoidal method estimates y., 


ee from the formula 


BC oe . 

V4 4% = YG + ope (y; * Yaa: (27) 
Substituting for Ys and Youd from Eq. 25 and rearranging 
gives 

Viel | 2-- dat 


y, 2+ Adt’ (28) 


Series exyansion of the right-hand side provides an error 


bound (per step): 
3 
(AAt)~/12. 


From the point of view of the size of AAt this method has no 


stability limit, but has slow attenuation with alteration 
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TABLE { 
ATTENUATION FACTOR COMPARISON 


Fourth Order Runge-Kutta Alvorithm 


At Yao /My exp (-AAt) 
(Runge-Kutta) (Exact) 
.0001 9599 8999 
.001 .9990 9950 
.01 .9901 9901 
Pai .9048 9048 
<2 .8187 .8187 
“ . 6068 .6065 
1.0 .3750 - 3679 
2.0 ¥5995 .1353 
25 6484 .0821 
3.0 1.3750 0498 
4.0 5.0000 .0183 
8.0 110.3333 0003 
10.0 291.0000 .0000 
20.0 §514.3333 0000 
50.0 240784.3333 0000 
100.0 4004901.0000 .0000 
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in sign for large AAt. Table II shows this behavior. 

Since a wide usable range of AAt is essential and 
the stability of trapezoidal integration is guaranteed, 
this method is chosen for the present program. 


Applying the trapezoidal algorithm to Eq. 13 yields 


hae ge ae) (29) 
where 
A=c+ Sty’, 
ip. At yt 
G c ae a 


and the superscripts denote evaluation at discrete time 
intervals At. If mis the order of the capacitance and 
admittance matrices, € and Y, then once a certain step size 
At is chosen, it requires m?/3 operations to perform the 
needed tx cngular decomposition of A. Thus, for large m, a 
change of step size At becomes costly from the point of view 
of computer time. Accordingly, in the present program only 
one time step size is used throughout each problem. 

Also, for assuring sufficiently rapid attenuation of the 
components corresponding to the large AAt, the following 
correction is utilized. 

“. Irons' Correction 

Irons proposed a scheme [4] for augmenting the atten- 
uation of the contributions of those eigenvectors for which 


Adt is large. Define 
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TABLE II 
ATTENUATION FACTOR COMPARISON 


Trapezoidal Integration 


re Yae/%G eg tht 
(Trapezoidal) (Exact) 

.0001 ~9999 .9999 
.001 -9990 .9990 
.01 9901 -9901 
ey -9048 .9048 
«2 - 8182 8187 
eS 7391 7408 
25 - 6000 .6065 
1.0 1oe55 . 3679 
2.0 .0000 _ 1353 
2.5 -.1111 0821 
3.0 -.2000 0498 
4.0 -.3333 0183 
8.0 - 6000 .0003 
10.0 - 6667 -0000 
20.0 -.8182 0000 
50.0 ~.9231 -0000 
100.0 - 9608 .0000 
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; 
} 


1 
} 
i 
H 
i 
j 
: 
4 


y= eeo Yz-3 7-5 vet -25 Yiyq> (30) 


where 3 and Yi, are obtained from Yz-4 by trapezoidal 
integration. 

In the program presented in Appendix B, Eq. 30 is used 
after every 10 steps of time integration. Table III shows 


the resulting modifications. 


{ 
E. ESTIMATION OF EXTREME EIGENVALUES 


Analytic results for one-dimensional heat conduction 


give, for an eigenvalue, 


A" toe iT? (31) 


where d is the distance between points of extreme tempera- 
ture and zero temperature. 
If we use this to estimate the smallest X in cylindrical 
coordinates, two modifications are recommended. 
1. Assume that the point of zero temperature is in 
the fluid at a distance from the wall equal to k/h, 
where h is the surface heat transfer coefficient. 
2. If there are two approximately orthogonal paths for 
heat flow from the (single) maximum temperature point, 


then replace 1/4? in the above formula by 


Bi, got dn (32) 
a2 


‘min. max. 


1 = 
a? 


For estimating the largest \, the surface heat transfer 


coefficient has no significant effect. We may continue to 
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TABLE III 


EFFECTS OF USING IRONS' CORRECTION 
AFTER 10 STEPS OF INTEGRATION 


Te exp (-10AA7) Yi0/Yo Ye aati 
Exact Trapezoidal Corrected 
0.00010 0.99900 0.99900 : 0.99900 
0.00020 0.99800 0.99800 0.99800 
0.00100 ) 0.99005 0.99005 0.99005 
0.01090 0.90484 0.90484 0.90486 
0.10000 0.36788 0.36757 0.36849 
0.20000 0.13534 0.13443 0.13579 
0.30000 0.04979 0.04866 0.04978 
0.50000 0.00674 0.00605 0.00645 
1.00000 0.00005 : 0.00002 0.00002 
2.00000 0.00000 0.00000 | 0.00000 
4.00000 0.60000 0.00002 — -0.00001 
8.00000 0.00000 0.00605 -0.00040 
10.00000 0.00000 0.01734 -0.00072 
20.00000 0.00000 0.13443 -0.00136 
50.00000 0.00000 0.44914 -0.00072 
100.00000 0.00000 0.67028 -0.00027 


use the same formula for 4 but now consider only the 
d 


smallest element and take 


a. = Length of smallest side 


min 4 > 
(32") 
d =. length of largest side 


max 4 


A comparison of estimates based on Eq. 31, 32, 32' with 
the exact solution of the eigenvalue problem has been carried 
out for several examples. Based on these comparisons, it is 
believed that these estimates are sufficiently accurate for 
choosing a time step and estimating the time of occurrence 


of the extreme stresses. 
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III. ONE-DIMENSIONAL HEAT CONDUCTION 


For comparison of numerical (time) integration methods, 
studies of one-dimensional heat conduction were made. In 
this section numerical results for trapezoidal time integra- 
tion using Irons' correction are compared with the exact 
transient temperature solution. 

Consider a flat slab of thickness L with conductivity k, 
density p, specific heat c, zero initial temperature, one 
face insulated and the other in contact with fluid at tem- 
perature Tr (Fig. 1). The surface heat transfer is h. The 


exact transient heat conduction solution is available [5] as 


aes. Sinu.L Cosp.x -y2 xk, 
T=Te]1-22 oypet tos | 
i=1 "i Bi ar 


Fluid 
Temperature T¢ 


— — 


—_ — 


Figure 1. Slab with One Face Insulated and Another 
in Contact with Fluid. 
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where k/pc is the diffusivity and the u; are the solutions 


of the transcendental equation 


= oe nine ect ermine remetrtetie endl Dinah wench 


hh, 
Tan wL = a TL e (34) 


For the finite element comparisons we subdivide the 
distance L into m one-dimensional 3-noded elements and use 
the corresponding isoparametric shape functions. (The work- 
ing equations, shape functions and the element capacitance 


and admittance matrices are given in Appendix A.) 


‘ ves 
+ Nmap boey Sosa alratd staan taut enlialonen ieee ikea dabei esiahi isda hat bala 


In the course of this investigation separate computer 
programs were written to evaluate nodal transient tempera- 
tures using the following methods: 

(a) exact transient temperature solution, Eq. 33; 


(b) exact time solution with spatial discretization 
? (section II.C, Eq. 24); 


Pho POEL A ay EMOTE A BT A Te aT? ¥ 


(c) Runge-Kutta time integration (section II.D, 
, part a); 


(d) trapezoidal time integration (section II.D, 
part b); 


(e) trapezoidal time integration with Irons' 
correction (section II.D, part c). 


For an initial step change of fluid temperature from 
zero to 1, transient temperatures have been found. For these 


comparisons the parameters (in consistent units) were taken 


2A SRR SRE RRR SEAN RRERE TNE 


to be: 
L = 8, p = 25, k= 8, c= 5 andhe= 5 


‘ The distance L was subdivided into m.= 3 elements. For the 
present purpose, comparison is confined to the exact solu- 


tion (item (a)), and the finally adopted system (item (e)). 


| 
| 
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In Table IV, temperatures at the two faces (x = 0, x = 8) 
are compared for various times. The trapezoidal integration 
has been performed using the constant time increment unity 
and the Irons' correction is applied after every 10 incre- 
ments. 

It is believed that Table IV demonstrates that the 
numerical integration method gives adequate accuracy for 


engineering applications. 


TABLE IV 


COMPARISON OF ONE-DIMENSIONAL TRANSIENT TEMPERATURES 


Exact 


Trapezoidal 
with Irons' 


Exact 


Trapezoidal 
with Irons' 


v 


Trons' correction is used after every 10 steps of trapezoidal 
integration. 
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IV. STRESS PROBLEM 


aR 
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For bodies of revolution deformed symmetrically under 
axisymmetric loading, the stress distribution is two- 


dimensional. Since deformation is symmetric about the axis 


Bence. Ala weg ne EN 
‘ 


Ts 
att 


abba 


BSR re 


of revolution, cylindrical coordinates (R,Z,0) are used. 

It follows that the stress components are independent of the 
angle 6 and all derivatives with respect to 6 are zero. 

Also the components of shearing stress Tre and Tz vanish 


on account of the symmetry. But since any radial displace- 


ee Ae 5 OREN COI 


ment induces a strain Eg in the circumferential direction, 
this non-zero component of strain and the three in-plane 
components (e7s Ep» Ypz)> complete the state cf strain at a 
point in any axisymmetric situation. Hence the state of 


, Stress for an axisymmetric body under axisymmetric loading 


ot iy aR ee, TE 
e 


is given by 


= sy 
G = <Oy Op Gg Thy? . (35) 


2 gt OR reas ae ROY 


In this chapter the stiffness matrix of an axisymmetric 


peed 


body and the thermal, pressure, and centrifugal load vectors 


are formulated and, finally, evaluation of stresses at a 


Sa are rae 


point is discussed. The treatment closely follows that of 


ARS 


Zienkiewicz [1] and this reference should be consulted for 


further details. 
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A. STIFFNESS MATRIX 

The elements used are bodies of revolution (about the 
Z axis). For analysis it is sufficient to describe the 
cross-section in the R,Z plane. In Fig. 2 such an element 


and the local §,n coordinates are shown. 


Figure 2. Quadrilateral Element Representation. 


If u and w are the displacement components at a point in the 
directions of R and Z respectively, then these displacement 
components may be defined in terms of the nodal displacements 


by the appropriate isoparametric shape functions as 
8 8 
Me oS Nelle 4 owe YE Nes -g (36) 


where Nj, a function of € and n , is the shape function for 
element node i and u;.W;, are the nodal displacement compo- 
nents. The strain-displacement relations [6] can now be used 


to obtain the components of the strain vector. Thus 
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€¢ = 


pas 


IIs 


& ; (37) 


where B is the standard rectangular strain-displacement 
matrix of any finite element formulation, a function of 
the local coordinates € and n, and § is the vector of nodal 
displacement. (See Appendix A, part 3, where the applicable 
formulas and useful equations are summarized). 

If the elasticity matrix for an isotropic material is 


D, then the stress vector o at a point is given by 
oe (38) 


Now, by evaluation of the total strain energy in the element, 


the element stiffness matrix can readily be obtained as 


DBa, (39) 


where the integration extends over the volume of the 
. element. 
In the present program the upper triangle of each ele- 
ment stiffness matrix is evaluated by numerical integration 
using four Gauss points within the range of &€ and of n [1]. 
The element contribution is immediately placed in the system 


stiffness matrix, which is stored in banded form. 


B. THERMAL LOAD VECTOR 
If we denote T as the difference between local tempera- 


ture and reference temperature, then the thermal strain Ey 


is given as 


= UoT, (40) 
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T 


U = <1 11 «0> 


and a is the coefficient of thermal expansion. 
The thermal load vector Le is given by Zienkiewicz [1] 
as 


an BI 


D CM dv. (41) 


From the point of view of the numerical evaluation it is 
interesting to note, however, thet the product D oe in 


Eq. 41 will reduce to 


De® = ae (42) 


where E is the modulus of elasticity and v is Poisson's 


ratio. Thus 


e | 
fp aay 


ua, (43) 


where T = ‘i N.T. and n is the number of nodes in each 
element. : 

In the attached computer program in Appendix B the advan- 
tage of the simplicity of Eq. 42 has been utilized. Also, 
since B has some zero components, in the process of multi- 
plication of BI U, simply the addition of the appropriate 
non-zero components of each column of the Bp has been per- 
formed. Finally, Eq. 43 has been integrated with four Gauss 


points. 
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C. PRESSURE LOAD VECTOR 

Consider a quadrilateral element as in Fig. 3 on the 
boundary of the axisymmetric cross-section where constant 
pressure P is applied. The infinitesimal force dF due to 
the normal pressure acting on the inner infinitesimal cir- 


cumferential surface dS is 
dF = PdS = 27RP dg , (44) 


Where d&’ is the infinitesimal length along the side of 


quadrilateral. 


Rv 


Z,w 


Figure 3. Boundary Element Under Pressure. 


Let Fr and Fy be the components of the pressure force in the 


R and Z directions respectively, then 


dF, = dF Cos 6, dF, = -dF Sin 6, (45) 
where Sin @ = g and Cos 6 = a . Therefore 
dF, = (27 RP)dz , dF, = -(27 RP)dR. (46) 
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Since the total work W done by the normal force F is equal 


to the sum of the work done in R and Z directions, 


W 


Su dFp + fw dF, 


2uP (SR u dZ - SR w dR). (47) 


Now, by using the appropriate shape functions, each compo- 
nent of Eq. 47 may be defined in terms of the nodal values. 


Since 
e T e 
W= (6°) Fp, (48) 


where ES is the element pressure load vector 


vector contributed by node j is explicitly given as 


e i 
z FR 1 Lae 23 
Fo = = 2nP f (ENR) N; dé , (49) 
Jp Fy ot ON. 
j ea 2 
Pp 0g 


where N, in Eq. 49 is the appropriate shape function for 


node j. 


D. CENTRIFUGAL LOAD VECTOR 

Refer again to Fig. 2 and assume that the body is rotat- 
ing about the Z axis with aheular speed 2. Then the centri- 
fugal force per unit volume at a point distant R from the Z 
axis will be 02" R, where p is the density of the material. 


The work done in this case is 


W= s pn? Rua. (50) 
V 
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For constant pk” the corresponding element centrifugal load 
vector is readily obtained by evaluation of the components 


of the integral in Eq. 50 in terms of the nodal variables, 


1.€., 


. d&dn, (51) 


e Bt 2 
oe = 27p2" f fs (ZN; R;) ace J N, 


A Fak 


where det J is the determinant of the Jacobian coordinate 


transformation matrix (see Appendix A). 


E. STRUCTURAL BOUNDARY CONDITIONS 


AUTE 
GREER aR PEERS 


The structural boundary conditions implemented in the 


program are: 
(a) one or more nodes prevented from moving axially; 
(b) one or more nodes prevented from moving radially; 
(c) the right-hand end cross-sectional plane remains 


plane and the transmitted axial force is specified. 


Hereinafter this will be referred to as the plane- 
. end boundary condition. 


. tet sak ones 
_ So cag STR ES 
OF Re SNS ae Ae 
AL PERE AS 
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The computer program presented in Appendix B has the 
capability of handling any combination of the above men- 


tioned structural boundary conditions. For boundary condi- 


tions of the types (a) and (b), simply multiplying the 
corresponding diagonal component of the stiffness matrix by 


1° gives zero displacement (8] (for practical purposes). 


ee 
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For the boundary condition of type (c) both ends are initially 
fixed axially for all solutions. An additional solution is 

obtained for init axial displacement of one end. The axial 
force is evaluated for each solution. Superposition is per- 


formed by adding the displacement vectors for the given 
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loadings (thermal plus mechanical), plus an appropriate 
fraction of the vector found for unit axial displacement. 
H 
z 
This fraction is chosen so that the resultant axial load has 


the specified value. 


F. SYSTEM EQUATION SOLVER 

Once the desired structural boundary conditions are 
applied, then the problem is to find the nodal displacement 
vector 6, corresponding to a given number of load vectors. 


We have 


Il 


SF, (52) 


where K is the system stiffness matrix in banded form and F 
is a load vector. In the present computer program a single 
Cholesky decomposition is performed on K. Then, by a process 


of forward and back substitution, each load vector is 


replaced by the corresponding displacement vector. 


G. PRINCIPLE OF SUPERPOSITION 

Upon the evaluation of the displacement vectors due to 
the various types of loading, the principle of superposition 
can be applied on the displacement vectors. On each thermal 
displacement vector the displacement due to any other type 
of loading is superimposed and, as the result, the number of 
displacement vectors is reduced to the number of thermal load 


vectors. 


H. STRESS EVALUATION 
From the system displacement vector 6, tne displacement 


vector of each element 6° may be obtained easily. Then the 
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corresponding element strain vector e° at any point can be 


found from 


eS BS" (53) 
Finally, the corresponding element stress vector a” is 


obtained by 


go = D(s =.) « (54) 


Since normally the stresses on the inner and outer surfaces 
of axisymmetric bodies are desired, in the computer program 
presented in the Appendix B provision has been made to cal- 


culate the stresses at the two Gauss points corresponding 
1 


v3 
element (where n = +1). 


to — = ¢# on the inner and outer boundaries of each 
Upon the evaluation of the stresses at each point the 
mean stress and the octahedral shearing stress [7] are 
calculated. The program gives as output the extreme values 
of these stresses, the R and Z coordinates of the corres- 


ponding points, and the times of occurrence. 


V. ONE-DIMENSIONAL TRANSTENT STRESSES 


In this section a one-dimensional comparison of stresses 
is made between exact and finite element results. The tran- 
Sient temperature problem is the one previously described in 
Section III. 

If the slab edges are free to translate in the plane of 
the slab, but are prevented from rotating, the exact solu- 


tion for thermal stress [9] is 


Spey ee : 
oy o, = (Taye T3°% (55) 


where T is the local temperature (Eq. 33), and 


Tavg is the average temperature in the slab 


If we choose E = 2, a = .50, and v = 0 (all in consistent 
units), then the maximum stress obtained from the exact solu- 
tion is 


Ona -.477786 


and it occurs at x = 8,1 = 73. 
Using the finite element technique with trapezoidal time 
integration and Irons' correction every 10 steps, the maximum 


stress is found to be 
eee ; -.477778 


and it also occurs at x = 8, t = 73, as before. 


It is observed that the method chosen gives excellent 


results. 


VI. TEST PROBLEMS 


Program integrity and accuracy have been verified by 
solving a number of test problems. Since stresses, whose 
evaluation depends upon derivatives of displacements, are 
known to be less accurate than temperatures, comparisons 
with exact results are confined to stresses. Individual 


problems are described below. 


A. THICK CYLINDER 
Consider a thick cylinder with inside radius 30 inches 


and outside radius 50 inches and the following material 


properties. 
Modulus of elasticity ; E = 28.9 x 10° Psi 
Poisson's ratio v= .28 . 


Coefficient of thermal expansion o = 7.22 x 107° 1/F° 


Thermal conductivity k = 28. Kee 
Density p = 489, Lbm/ft> 
Specific heat c= .1lil ets 


An arbitrary length of 25 inches has been selected for the 
cylinder and it has been subdivided into two different 
element representations as in Figs. 4 and 5. The plane-end 
boundary condition with zero axial force is used. The 
stresses for various types of loading are compared with the 


corresponding exact analytic solutions as described below, 


Figure 4. Two Radial Elements Representation 
of Thick Cylinder. 


Figure 5. Five Radial Element Representation 
of Thick Cylinder. 


1. Thermal Loading 
For a linear variation of the temperature T = 20R 
the stresses obtained by the finite element niethod: for the 
above representations are compared with the exact analytic 
solution in Fig. 6. The Tp for this problem clearly is 
zero and the one obtained by the program was 10", The 
accuracy of the other results is clearly satisfactory. 
2. Pressure Loading 
A uniform pressure of 1000 psi. acts on the inner 
surface of the cylinder. Again, Tee is zero and the program 
gives 19°19, In Fig. 7 the other stresses induced by this 


uniform pressure are compared with the exact solutions. 


Here aiso the accuracy of the results, even with two radial 


elements, is adequate. 
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Fig.7 THICK CYLINDER UNDER INTERNAL PRESSURE, 
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3. Centrifugal Loading 


The speed of rotation has been assumed to be 500 
revolutions per minute. The stresses obtained by the pro- 
gram are compared with the analytic solution in Fig. 8. In 
this case also the results obtained by the program, even 
with only two radial elements, are very close to the exact 


solution. 


B. HOLLOW SPHERE 

We consider a hollow sphere with the same material 
properties as in the thick cylinder with inside spherical 
radius 30 inches and outside 50 inches. The loading is 
thermal with T = 20 ‘(spherical radius). Symmetry permits 
using only half of the sphere for the computer analysis. 
The elements representation is given in Fig. 9. Since the 
program gives stresses in cylindrical coordinates, these 
have been transformed to the spherical coordinates for 
comparison with the exact solution in Fig. 10. The accuracy 


of the results is noteworthy. 


C. THERMAL STRESSES IN NOZZLE 

This problem concerns thermal stresses near the inter- 
section of a cylindrical pipe and the spherical vessel. 
Fig. 11 gives the cross-section of the structure which is 


to be analysed. The material properties are: 


Modulus of elasticity E = 29.3 x 10° Psi 
Poisson's ratio | v= .30 - 
Thermal expansion coefficient a = 7.6 x 10°° 1/F° 


25 —~ EXACT 
B 2 ELEMENTS. 
2. o 5 ELEMENTS. 
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FIG, 8 
ROTATING CYLINDER 


Figure 9. Hollow Semi-sphere 32 Elements Representation. 
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Thermal Stresses in Hollow Sphere 
T = 20 (spherical radius). 
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Thermal conductivity k 10.25 Btu/ft.hr.°F 


530.5  Lbm/ft~: 


tt 


Density p 


: Specific heat c 128 Btu/Lbm°F 


The loading results from the thermal transient in the fluid 
contained in the nozzle. This fluid is in contact with the 
structure on surface "1" (Fig. 11) and has the entry temper~- 
ature time variation as in Fig. 12. The fluid in the sphere, 
which is in contact with the structure on surface "2" (Fig. 
11), has the constant temperature 478°F. At t = 0 structure 
has a uniform temperature of 478°F and is stress free. 
Exterior surface of the structure is susulated: The flow 
velocity past surface 1 is 8.5 ft/sec (inward) and 

there is no flow past surface 2. The surface heat 

transfer coefficients are 1393 and 2910 Btu/hr. £t2°F for 
surfaces "1" and "2" respectively. 

For the structural boundary condition it is assumed 
that the nodes on the left end cross-section of Fig. 13 are 
prevented from any axial displacement. 

The maximum thermal stresses obtained by the program 


occur in element 14 of Fig. 13 as follows: 


Maximum mean stress = 18.81 ksi. at time = 4 seconds; 


Maximum octahedral shearing stress = 11.5 ksi. at time 
= 18 seconds. 


‘ These results appear to be reasonable, but no suitable com- 


parison solution is available. 


at 


478 


ee eee 
O 2 SEC. 
VARIATION OF ENTRY FLUID TEMPERATURE 


IN NOZZLE, (SURFACE ‘!".) 


478 


O A re rrr fe 
SEC. 
FLUID PAST SURFACE “2” 


Fig. 12, Fluid Temperature - Time Histories. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

A computer program has been developed for the solution 
of axisymmetric transient heat conduction and thermal stress 
problems. The system will accommodate a wide variety of 
geometric arrangements, thermal and structural boundary con- 
ditions, and mechanizal loadings. 

Although double-precision arithmetic is employed through- 
out, efficient algorithms for the manipulation and storage 
of large symmetric banded matrices allow in-core solutions 
with modest time requirements. 

The quadratic isoparametric elements used allow accurate 
representation of curvilinear boundaries and the stress 
field. Examples presented show that a small number of ele- 
ments is generally sufficient to determine stresses with 


good precision. 


B. RECOMMENDATIONS 

Incorporation of several additional features would 
Significantly increase the program capabilities. The follow- 
ing extensions are recommended. 

1. Material thermal and elastic properties have been 
assumed constant within each element. Since such properties 
are generally temperature dependent,.provisions should be 
made for periodically "updating" them during both temperature 


and stress solutions. 


2. The surface heat transfer coefficient, taken as 
constant in the program, is a function of temperature and 
flow velocity. Provision should be made to include thesé 
effects. 

3. The program presently starts every temperature 
solution with constant initial solid temperature. Provision 
for an externally specified initial temperature vector should 
be included. 

4. The large quantity of temperature and stress results 
generated by the program is currently presented as digital 
printout. Graphical output in the form of two-dimensional 
contour plots of temperatures and stresses should be 


provided. 
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APPENDIX A 
APPLICABLE FORMULAS AND EQUATIONS 


PART 1 


(a): Two dimensional 8-noded (parabolic) 
isoparametric shape functions. 


Corner nodes: 
1 
Ni 7 1+) + nC, tn, - YD 
Mid nodes: 


GQ - 64) +n) 


= 
SS) oe 


(+ & C1 - 07) 


3 
1" 
oO 
A 
W 
[SS 


where 


uty 
It 


b&3> n =nn, 


(b): Temperature at a point in terms of the 
nodal temperatures. 


(c): Coordinates at a point in terms of the 
nodal coordinates. 


(d): The Jacobian coordinate transformation 
matrix. 
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(e): Element of the finite element capacitance 
matrix. 


e 7 . 


(£): Element of the finite element admittance matrix. 


e | = os 
Yaz 7k S N,-WN, av 


(g): Elemental volume. 
dV = 2nR det J d& dn 


PART 2 


(a): One-dimensional 3-noded (parabolic) iso- 
parametric shape functions. 


End nodes N = . bo (1 + a) 
Mid node Nv = (1 - 5 
where, again, ae = b5 


(b): One-dimensional capacitance and admittance 


matrices. 
2 - 4 7 -8 1 h 0 0 
e _ pck 5 Saar Sa : 
Cc 30 2 16 res = 35 8 16 8} + 0 0 0 
- 2 1 -8 7 0 0 O 


where 2 is the length of the element. 


(c): The element v vector for the one-dimensional case. 


v° = <hT¢ 0 o>! 


PART 3 


(a): Strain-displacement relations. 
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where N are the same isoparametric shape 
function as in Part 1, (a). 


(b): The B matrix. 
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(c): The element displacement vector. 


e _ 
oF SU 


v Vv 

1 1l-v 1-v 0 

“ E(1-vj v 7 Vv 0 

= (itv) (i-2v) |iI-v l-v 
v v 

1l-v l-v 1 0 
0 0 1-2v 
2(1-v) 
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(d): Elasticity matrix D for an isotropic material. 
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APPENDIX C 
USER'S MANUAL 


In this section the procedure for the use of the present 
computer program is described. It is intended that a person 
with minimum familiarity with the details of the finite 
element method and computer programming be able to use this 
program. 

The steps below are in order and the user is advised to 
follow them carefully. 

For finding the thermal stresses in an axisymmetric body 
under axisymmetric loading, the user has the choice of using 
either the Metric or British system of units. The program 
will handle both systems and the necessary conversions are 
made automatically within the program. However, the units 
used in each system must be consistent and they should be as 
listed below in Table V, 

Now for preparation of the data input for a given axi- 
symmetric geometry with prescribed thermal and structural 
boundary conditions we go through the details with a simple *! 


example. 


Step 1: 
Draw the longitudinal cross-section of the body to scale. 
Identify the cylindrical coordinates R and Z, with the 


rigin of the Z axis on the most left-hand point of the 


TABLE V 
UNITS FOR INPUT DATA 


Note: an input card specifies whether British or Metric 
data is being used. 


Quantity British Unit Syst. Metric Syst. 
Coordinates in. cm 
Time sec. - sec. 
Temperature oF 2c 
Velocity ft./sec. m/sec. 
Axial force” lbf kg 
Pressure 1b£/in? kg/cm 
Density 1bm/ ft gm/cm? 
Specific heat Btu/1lbm°F cal/kg°C 
Mod. of elasticity 1b£/in? kg/mm? 
Coef. of thermal exp. L/°F 1/°c 
Thermal conductivity Btu/hr. ft. °F cal/sec.cm. °C 
Heat transfer coef. Btu/hr. ft. 2°F cal/sec.cm. “°C 
Load vector 1bf kg 
Rotational speed Revolutions/min, Revolutions/min. 


¥ 


lbf is pound force and lbm is pound mass. 
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body. The entire cross section must lie in the first quad- 


rant of the coordinate plane. 


R 
13, 
9 
5. 
° 5. 10. 15 Z 
Figure 14. Longitudinal Cross Section. 
Step 2: 


Subdivide the cross section into a maximum of 40 quadri- 
lateral elements. This subdivision is extremely important 
and we should provide smaller elements for the parts of the 
cross section where we expect the largest stresses or tem- 
perature gradient. If the body is made from several mater- 
ials, elements should be chosen so that each element contains 
only one material. Any two adjacent elements must share one 
complete side of the quadrilateral. Number the elements, 


starting from 1, in any arbitrary manner (see Fig. 15). 


Step 3: 


Since eight-noded elements are used in the program, iden- 
tify these nodal points around the boundary of each element. 


Four nodes will be at the corners and the other four will: be 
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ne ny mem nent 


at the mid-points of the sides. Number sequentially these 
nodes starting from 1 and increasing in the direction where 
the number of elements is the least. (fhe numbers thus 
assigned are called global node numbers.) For clarity of 
this step assume a cross-section as in Figure 15 such that 
the maximum number of elements in one direction is less than 
the maximum number of elements in other direction. (In 

Fig. 15 we have maxima of 2 and 3 elements in R and Z direc- 
tions respectively.) Thus we number the nodes beginning in 
the R direction. See Fig. 16. This method will give the 
minimum band width of the stiffness matrix and will save 


computer execution time. 


Figure 15. Subdivision Inot Elements. 


Step 4; 


At this point we are ready to prepare the first data 


card. Input quantities are: 


"Figure 16. -Element and Node Numbering. 


NET The total number of elements 
NNT Total number of nodes 
‘ NCN Total number of corner nodes 


3 NQ@EN Total nunber of mid-side nodes not lying on the 
am . straight line joining the corner nodes (number 
| of "Off" nodes) 


NMAT Number of different materials 


NPROB Number of problems to be solved 


P 4 For the present we take NPR@B = 1. 


: Prepare a Single card that reads 
NET, NNT, NCN, N@FN, NMAT, NPRQ@B 


with the format (815). For our example, if all the elements 


are from the same material this data card reads: 


'4 e@3 i { j ‘40 ; ee 


ig « : . b 
im Bahn eh SRT be to te Ute tin Se Sh ue ta tae Eh NTO HT Clam TR Bi bec,” sek 


Ce ow meee se erae Ce ees v4 Be Cake AL EMCI, GOA SAM ACS OE Beate Sed) AC AA oe 
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On this card, or any succeeding one, if an input quantity 
(e.g., N@FN) is zero, the corresponding field may have a 0 


or be left blank, unless otherwise stated. 


Step 5: 

For each corner node a card should be punched which 
seas: The node number I, the R coordinate and the Z coor- 
dinate of the node. The format is (I5, 2F10.5). 

It is not necessary to sort these cards in the order of 
increasing or decreasing corner node numbers, they can be 
put together in any arbitrary order. A typical card is 
illustrated in step 7. There must be as many punched cards 


as the number of corner nodes (NCN). 


Step 6: 

In this step we read — the connectivity array, i.e., 
for each element we prepare one card which gives the global 
node numbers and the material identification number. 

For each element, start from the lower left-hand node 
and move in the counterclockwise direction within that 
element and punch the global node numbers in order. The 
format is (915) where the last I5 is the material identifi- 
cation number. 

It is very important that the cards prepared for this 
step be put together in order of elements, i.e., the first 
card for element 1, the second card for element 2, etc. 

For ease of sorting the connectivity cards the element 
number may be punched after column 50. As an example, for 


element 2 of Fig. 16 the connectivity card should read: 
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where 1 identifies the type of material in this element and 


2 in column 55 stands for the element number. 


Step 7: 

If each side of every element is straight, NQGFN = 0 and 
this step is omitted. 

For each mid-side node that is on a curved element edge 
a card is prepared with format (15,2F10.5) to read the global 
node number and the R and Z coordinates. Here, as in step 5, 
these cards may be put together in any arbitrary order. For 
the case of Fig.16, there will be only one card to be punched 


and it would read: 
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To assure in all cases that the double precision capability 
of the program is not compromised by less precise floating 
point inputs, it is recommended that all such inputs be made 
in D-FORMAT. The F-FORMAT instructions merely allocate 


input card fields. 
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Step 8: 


The properties of the different materials are specified 
in this step. For each material a card must be prepared 
that gives the modulus of elasticity (E), coefficient of 
thermal expansion (AL), Poisson's ratio (P@I), thermal con- 
' ductivity (TK), density (DENS) and specific heat (SHT). 


I.e., read: 
E, AL, P@I, TK, DENS, SHT 


with the format (2D12.4, 4F8.3). The first card will be for 
material number 1, the second card for material number 2, 


etc. 


Step 9: 
For this step a single card must be punched, starting in 
column one, which reads the word BRITISH or METRIC in accor- 
dance with the system of units used. 
This completes the input of geometric and maierial prop- 


erty data. 


Step 10: 


For a transient temperature problem only (no stress cal- 
culations) leave a blank card for this step and proceed 
directly to step 16. 

For a stress problem or thermal stress problem, in this 
step we specify the type of problem and the structural 
boundary conditions. 

The following quantities, when pertinent, are to be 


specified. 
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@MEGA The speed of rotation, about the Z axis, in 
revolutions per minute 


PRES The external pressure applied to a boundary 
segment 
TEXT The number of known temperature vectors for 


which evaluation of thermal stresses is desired 


L@AD The indication for an additional known load 
vector. If there is one, L@AD = 1, otherwise 
LGAD = 0 


NNLT Total number of nodes fixed in the longitudinal 
direction 


NNRT Total number of nodes fixed in the radial 
, direction ; 


IPLANE An indication for the plane-end boundary condi- 
tion. If desired IPLANE = 1 


If any of these items is not applicable, the correspond- 
ing field is left blank. 
Now we have to punch a single card for this step which 


reads: 
(MEGA, PRES, IEXT, LOAD, NNLT, NNRT, IPLANE 
with format (2F10.4, 515). 


Step li: 


If there is no pressure loading (PRES=0), omit this 
step. For non-zero pressure, the total number of nodes on 
pressure loading boundary segment (NPNT) and the global 
node numbers of these (pressure) nodes (NPN(I)) must be 
specified here. Only one segment is permitted. 


We read in: 
NPNT, (NPN(I),1=1,NPNT) 


with the format (1215). 
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In our example, Fig. 16, if there is pressure inside the 


body, then the data card for this step will be: 


roa pales ia Nem aS as 5 apig Na: . “ arch aa top" OTIS ER BN TIE I to og OF 
ee Debi S BARROS AGL eS AER EEE 6 ome 14 aE NSE LOPE ATS e TIN Pe ce Pee oe 


OCOOOONDDODDNODNDNDOOHONFODNNN ON OO OCR HNNR RH ECOONCO ECOL 
9 NOD 1213 1615 161718 19 20.20 2223242 25 2 28 25 29.152 33395 1530 938 £9 2 045A A GE 525954 559657 SESE EO ENG 


er ws Cseevsveecer ger travresreesegars 


CJ 
LJ 
was 
3 
J 
D2 
ba —a 
. J 


- Step 12: 

If there is no additional load vector (LPAD = 0), omit 
this step. 

For LGAD = 1 read in the additional load vector with the 
format (6F10.4). The components of the load vector are 
arranged in the order: R component at node 1, Z component 
at node 1, R component at node 2, etc. JI.e., READ 
(F(I) ,I=1,NDF), where NDF (the number of degrees of freedom) 


is equal to two times the total number of nodes (2*NNT). 


Step 13: 

In any stress problem, or thermal stress problem, there 
must be at least one node constrained against longitudinal 
motion, i.e., NNLT>0O, so we must specify here which nodes 
are to be constrained against longitudinal motion. These 
global node numbers NNL(I) are read in with the format 


(1015), i.e., 
READ: (NN(I),I=1,NNLT) 


In our example, if global nodes 1 and 21 are fixed in the Z 


direction, we have a card that reads: 
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Step 14: 

If there are no nodes fixed in the radial direction 
(NNRT = 0), omit this step. 

For NNRT>0O we must read in the global numbers (NNR(I)) 
of those nodes which are to be fixed against radial motion. 


The format is (1015), i.e., 
READ: (NNR(I),I=1,NNRT) 


If node 1 and 21 of our example are also fixed in the R 
direction, then the card for this step would read exactly 


as the one in Step 13. 


Step 15: 
If there is no plane-end boundary condition (IPLANE = 0), 


omit this step. 

For the case of end-planes-remain-plane boundary condi- 
tion, i.e., (IPLANE = 1), we have to specify the total number 
of nodes of the right-hand end (NNRE) and the global node 
numbers of this end (NNR(I)). 


We read in: 
NNRE, (NNR(I),I=1,NNRE) 


with the format (1015). 
For the case of Fig. 16, if it is desired to have end- 
planes-remain-plane boundary condition then the data for 


this step will be: 
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Step 16: 


In this step we begin the specification of thermal 
boundary conditions. Additional details are prescribed in 
Steps 19 through 24. 

Thermal boundary conditions are imposed along discrete 
segments of the boundary. Each such segment must begin and 
end at a corner node of a boundary element. We consider 
separately the imposition of the various kinds of thermal 
boundary conditions. 

(a) Convection 

Since the fluid temperature for convection is deter- 
mined from a prescribed temperature-time history at entry 
section and a specified flow velocity (see Appendix D, 
Section 5), it is necessary that the terms “inside," i.e., 
adjacent to the symmetry axis, and "outside" have their 
ordinary meanings. Thus on the inside surface the convec- 
tion boundary condition may be applied to a single (con- 
tinuous) segment. A similar prescription may be employed 
for the outside portion. The entry section used for flow 
calculations is at the upstream end of corresponding segment. 

(b) Constant temperature 

To specify constant temperature on a portion of the 
boundary, the designator "inside" or “outside" may be used. 


Such a constant temperature portion may consist of several 


discrete segments. If two different constant temperature 
portions are prescribed, one may be designated "inside" and 
the other "outside." 


(c) Combinations of convection and 
constant temperature 


The convection and constant temperature conditions 
can be used together, but the portion designated "inside" 
must have only one of these conditions prescribed and the 
same restriction applies to the "outside." 

(d) Insulated 

All portions of the boundary not included in the 
segments specifically identified as "inside" and "outside" 
are considered insulated. 

The following items, when pertinent, are to be 


given as specified below: 


TINIT The constant initial body temperature 
Q2=0. For insulated boundary condition outside 


Q2>0. For convection or constant temperature boundary 
condition outside 


Q3=0. For insulated boundary condition inside 


Q3>0. For convection or constant temperature 
boundary condition inside 


Q4 <0. If solving stress problem only 
Q4>0. If solving temperature problem only 
Q4=0. If solving thermal stress problem 
AXIALF The axial force in the Z direction. As usual, 
+ (plus) for tension and - (minus) for compres- 
sion. ‘ 
So, we read: 


TINIT, Q2, Q3, Q4, AXIALF 
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with the format (5F10.4). 

In case of Fig. 16, if the initial solid temperature is 
60° and we have insulated outside and convection boundary 
condition inside wich 120 kg axial compressive force, then 


the card for this step will be: 
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where the blanks ieft in the columns 31 through 40 indicate 


that we want to solve a thermal stress problem. 


Step_17: 

No action is taken in this step - we merely choose between 
proceeding to step 18 or jumping to step 26. 

If there are no known temperature vectors for which 
evaluation of thermal stresses is desired, proceed to 
step 18. 

For IEXT>0, i.e., when some known temperature vectors 
are to be entered for thermal stress analysis alone or 
combined with some other loadings, proceed directly to 


step 26. 


If (C4<0); i.e., we are to solve only a stress problem, 


proceed directly to step 27. 
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Step 19: 


If there is an insulated boundary condition outside 
(Q2 = 0}, omit the following steps and start from step 22. 
For the case of a convection or constant temperature 
boundary condition outside (Q2>0.), we have to specify 
the outside heat transfer coefficient (HTC@} (for constant 
temperature HTC@ = i0°° ;}, the constant outside initial 
fluid temperature (TEMP®) which may be equal to the solid's 
initial temperature, the total number of nodes in contact 
with fluid outside (NCFPT), and the number of temperature 


ramps for outside flow (NRAMP@). We read: 
HTC@, TEMPO, NCKAT, NRAMPG 


with the format (D16.4,Fi0.4,215). See example in step 22. 


Step 20: 

Here we read in the gloval node numbers of the nodes 
in contact with outside fiuid (NCF%(I))}. The sequence of 
these node numbers must be in the direction of flow velocity. 


Read: 
(NCFG(1) , I=1,NCFOT) 
with the format (1215) (see example in step 23). 
Step 21: 
For each ramp of outside flow we read in the time when 
the ramp starts BDRY@(I,i), the initial temperature of the 


ramp BDRY#(1,2), the f*nal temperature of the ramp (specify 


only if it differs from the initial temperature of the next 


ramp} BDRY#(1I,3) and the velocity of the fluid for this 
ramp BDRY@(I,4), with the format (4710.4) so we read: 


((BDRY@(1,J} ,J=1,4) ,I=1,NRAMPQ) . 


The number of cards prepared at this step is equal to the 


“number of ramys. See example in step 24, 


Step 22: 

If there is an insulated boundary condition inside 
(Q3+0),. proceed directly to step 25. 

For the case of convection or constant temperature 
boundary condition inside, Q3>0. We have to specify the 
inside heat transfer coefficient (HTCI) (for constant 
inside temperature HTCI=102°) , the constant inside initial 
fluid temperature (TEMPI), which may be equal to the solid's 
initial temperature, total number of nodes in contact with 
fluid inside (NCFIT) and the number of ramps of the inside 


flcw (NRAMPI). We read: 
HTCI, TEMPI, NCFIT, NRAMPI 


with the format (016.4, F1l0.4,215). 

For our example assume th2 time variation of entry temper- 

ature for the inside flow to be as in Fiy. 17. The heat 

transfer coefficient is 175.0 with initial inside fluid 

temperature 70.0, The card for this step is: 
" "(75.00 rele ? 3 
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Figure 17. Entry Temperature-Time Variation. 


Further details of this temperature-time history are speci- 


fied in step 24, 


Step 23: 


Here we read in the global numbers of those nodes which 
are in contact with inside fluid, (NCFI(I)) with the format 
(1215). It is important to read these node numbers in the 
direction of flow velocity. 

For example, if the inside flow of Fig. 16 is from left 


to right, then for this step the card reads: 
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Step 24: 


In this step data for inside flow are given. For each 
ramp I we read: the time when the ramp starts BDRYI(I,1), 
the initial temperature of ramp I BDRYI(I,2), the final 
temperature of ramp I (specify only if it differs from the 
initial temperature of the next ramp) BDRYI(I,3), and the 


velocity of the fluid for this ramp BDRYI(I,4). We read in: 
((BDRYI (I,J) ,J=1.4} ,I=1,NRAMPI) 


with format (4F10.4). 
For example of Fig. 17, let the flow velocity for the 
first ramp be 15 and that for the second and third segment 


be 20. Three cards are needed: 
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Step 25: 


In this step we must specify the following items: 


DTI The time increment (step size) of trapezoidal 
integration 


TIME1 The time when the first calculated temperature 


vector is to be stored for thermal stress 
evaluation 
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EVERY The interval of time between the successive 
temperature vectors which are to be stored 


IVEC Total number of temperature vectors to be 
stored for thermal stress evaluation 


INTP The indication for printing the calculated 
temperature vectors 


If INTP = 0 there will be no temperature prints 


If INTP = 1 the program will print the tem- 
perature after every step of time integration 


a If INTP = 5 the program will print the tem- 


a perature after every 5 steps of time inte- 
a : : grations, etc. 


At this step we prepare a single card that reads: 


DTI, TIME], EVERY, IVEC, INTP 


with format (3F10.4,2I5). 


\ 
For example, if the step size of trapezoidal time inte- 


gration is 2 sec. and we desire to evaluate thermal stresses 
for 10 different temperature vectors, each 40 seconds apart, 
and starting with the first thermal stress calculation at 


time equal to 60, then the card for this step is: 
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where 15 in columns 34 and 35 indicates that the temperatures 
will be printed every 15 steps of time integration. The 
total time of integration for, this example is 420 seconds, 


since 


TIME] + EVERY * (IVEC-1) = 420. 


Step 26: 


If there are no known temperature vectors (IEXT = 0) 
for which thermal stress evaluation is desired, omit this 
step. 

For IEXT>0O we must specify the nodal temperatures 
which are to be stored for thermal stress evaluation. We 


read: 
((STOR (I,J) , J=1,NNT), I=1,1EXT) 


with the format (6F10.4). 

This means that we enter all the components of the first 
nodal temperature vector, followed by the components of 
the second nodal temperature vector and so on until JEXT 


such vectors have been entered. 


Step 27: 


The data cards for this problem are completed now. If 
no more problems are to be solved for the same geometry, 
omit this step. 

If another problem is to be solved for the same geometry 
and spatial discretization, increment by 1 the NPR@B in step 


4 and start the input of the new problem from step 10. 


Step 28: 

If there is no access to the IBM 360 computer of the 
Naval Postgraduate School, omit the following steps and 
Start from step 32. 

For the convenience of the user, the author has put 


a so-called CHECK program and the program presented in 


Appendix B in the computer system at N.P.S. (Naval Post- 
graduate School). It is advised, however, that the user 
check his input data deck with the CHECK program before 
attempting to solve the problem. 

For using the CHECK program prepare the following 


control cards. 


//XXXX0000 JMB (0000, 0000FT,XX00), 'NAME' , TIME=1 

//JOBLIB DD DSN=F0609.BAKH, DISP=S'iR, UNIT=2314 , VOL=SER=DUFFY 
//G® EXEC PGM=CHECK,REGION=100K 

//FTO6F001 DD SYSQUT=A, DCB=(RECEM=FBA, LRECL=133, BLKSIZE=3325), 
// UNIT=SYSOUT 

//FTOSFO01 DD # 


where the first card is the regular FORTRAN job card used 
at this institution. 
Prepare your deck as Fig. 18 and it may be read in from 


the hot card reader. 


Data deck 


Control cards of step 28 
Job card 


Figure 18. Set-up for Using CHECK Program. 


If the output of the CHECK program has any error messages 
or has not been run, then there are some mistakes in the 
data deck or control cards. 

If there are no error messages on the output of CHECK 
program, then study carefully the output and make sure it 
is the same problem that has been intended to be solved. 
Once the correctness of the input data has been insured, 


proceed to the next step. 


Step 29: 

AXITTS stands for Axisymmetric Transient Thermal Stress 
and this is the name given to the program presented in 
Appendix B when stored in the computer. For using that we 


must prepare the following control cards. 


//XXXX0000 JOB (0000,000FT,XX00) 'NAME' , TIME=10 

//JOBLIB DD DSN=F0609. BAKH, DISP=SHR,UNIT=2314 , V@L=SER=DUFFY 
//GG EXEC PGM=AXITTS, REGION=425K 

//FTOGF001 DD SYS@UT=A,DCB=(RECFM=FBA, LRECL=133, BLKSIZE=3325), 
// UNIT=SYSQUT , SPACE= (CYL, (6,1) ) 

//FTOSFO001 DD * 


where the first card is the regular F@RTRAN job card. It is 
advised to ask for 10 minutes time, i.e., TIME=10, since 
this would not affect the priority of the job within the 
class K jobs. 

Prepare the deck as in Fig. 19 (it may be read into the 


computer from the hot card reader) and submit a so-called 


service request card, available in the computer center, to 


the operator on duty. 


Data deck 


Control cards of step 29 
Regular FGRTRAN job card (TIME=10) 


Figure 19. Deck Set-up for Using AXITTS Program. 


Step 30: 


If the user wishes to have a listing of the program he 
May prepare the following control cards as in Fig. 20 since 


the program is also listed on the data cell for this purpose. 


Regular F@RTRAN job card 
//PRNT EXEC PGN=IEBPTCH 
//SYSPRINT DD DUNMY 
//SYSUT1 DD DSN=F0609. BAKH, DISP=(LD, UNIT=2321 , 
//  V@L=SER=CEL003, DCB= (RECFM=FB, LRECL=80, BLKSIZE=2000) 
//SYSUT2 DD SYSPUT=A, SPACE=(CYL,6) 
//SYSIN DD * | 
PRINT TYP@RG=PG,MAXFLDS=1,MAXNAME=1 
MEMBER * NAME=AXITTS 
RECORD FIELD=(80) 
[* 


Figure 20. Deck, Set-up for Obtaining a Listing 
of the Program AXITTS. 
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The set-up deck of Fig. 20 may be read in from the hot card 


reader to get a listing of the program. 


Step 31: 
If the user wishes to obtain a deck of the program AXITTS 


he may prepare a deck as in Fig. 21 and read it in from the 
hot card reader. He must notify the operator on duty to be 


prepared for punching almost two boxes of IBM cards. 


Regular FORTRAN job card 
//PUNCH EXEC PGM=IEBPTCH 
//SYSPRINT DD DUMMY 
//SYSUT1 DD DSN=F-0909.BAKH,DISP=@LD,UNIT=2321, 
// V@L=SER=CEL003, DCB= (RECFM=FB, LRECL=80, BLKSI ZE=2000) 
//SYSUT2 DD SYS#UT=B,SPACE= (CYL, 6) 
//SYSIN DD #* . 
PUNCH TYP@RG=PQ, MAXFLDS=1 , MAXNAME=1 
MEMBER NAME=AXITTS 
RECGRD FIELD=(80) 


/% 
Figure 2]. Deck Set-up for Obtaining a Punched 
Deck of the Program AXITTS 
Step 32: 


For using the program presented in Appendix B, if the 


user does not have access to the computer facility at N.P.S., 


he must punch a copy of the program.: (Good luck!) 


For the stcrage location 425K bytes are required and 
since the output may be longer than what usually is allowed, 
three additional cylinders are recommended. The execution 
time required depends on the size of the problem and the 
number of time integration steps; however, 10 minutes 
computer time would be sufficient for a problem of 140 nodes 
and 500 steps of time integration, with 20 temperature 


vectors stored for stress calculation. 
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APPENDIX D 
PROGRAMMING 


The computer program presented in Appendix B is formed 
from one main program and fifteen different subroutines. 
The communication between the main program and the sub- 
routines is handled through the common blocks. In order to 
minimize storage requirements, most of the storage location 
of the system stiffness matrix is used for temperature 
calculations. This is accomplished by use of EQUIVALENCE 
seavenents. The total storage requirement is 425K bytes. 

In this section the function of some parts of the pro- 
gram is discussed and the assumptions used are brought to 


attention. 


1. MAIN PROGRAM 
The main program is simply calling different subroutines 
when they are needed. The subroutines which are called in 


main program are; 


INPUT, PR@B, CANDY, FLOW, FGRMV, TEMPER, PLOT, STIFF, FORMF, 


CENT, PRESS, DISPL, AND STRESS. 


2. SUBROUTINE INPUT 
In this subroutine the data regarding the geometry of 
the body, subdivisions into elements, and the material 


properties are read in and printed out. Calculation of the 


mid-side coordinates based on the straight line is done 


* ta Fa, Ore ele ie 


here. This subroutine calls for subroutine PLYTRZ once 

to produce a graphical representation of the nodal points 
then the half-band-width of the system stiffness matrix is 
calculated from the connectivity array. Also, depending 
upon the choice of the system of units, the necessary con- 
versions in each system of units are accomplished and the 
reference temperature (70°F or 20°C) based on the system of 


units is selected here. 


3. SUBROUTINE PROB 

For each problem this sutroutine is called once by the 
main program to read in and print out the information 
regarding the nature ef the problem. Based upon the given 
information the type of the problem is distinguished here 
and for the thermal problems the length of time integration 


is calculated. 


4, SUBROUTINE CANDY 
Subroutine CANDY (C and Y) evaluates the capacitance 

matrix C and admittance matrix y" of Eq. 13, in banded form. 
For each problem this subroutine is called once. The func- 
tional flow chart of the subroutine CANDY is given in 

Fig. 22. Since the only non-zero elements of the symmetric 
matrix Y" (Eq. 11) are the diagonal, first and second super- 
diagonals, in order to conserve storage and reduce the num- 


ber of arithmetic operations, the non-zero elements are 


stored in three different vectors (DIAG, @FF1, @FF2). 
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Figure 22. Functional Flow Chart of CANDY. 
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Once the system capacitance and admittance matrices C 
+ ie . 5 
and Y are formei, then the matrices A and G (Eq. 29) are 


evaluated and replace C and Y respectively. 


S. SUBR@UTINE FLEW 

Ix there is any convection or constant temperature 
thermal boundary coadition, this subroutine is called from 
the main program for each step of time integraticn ia order 
to evaluate the temperature of the fluid nodes. The calcu- 
lation is based upon the constant fluid flux assumption 
(for beth inside and outside flow). Consider a section of 
an irregular cylindrical pipe as Fig. 23 and focus attention 


on eremenc itl on the inner boundary. 
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Figure 23. Fluid Node Representation for Inside Flow. 


If Fluid is moving from left to right we may number the R 
and Z coordinates cf the fluid nodes of the element itl as 
in Fig. 23, then the volume of fluid pumped into the pipe 


at time t is given by 
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Q;.. = ave, (56) 
in 


where a and v are, respéntively, the area and velocity at 
the entrance section. Also the volume of the pipe up to 


the mid-node of element itl is 


Q = V; + f mR” dZ. (57) 


where V; is the volume of the pipe between the entrance 
section and element i+l on the inner boundary. 


Now we can write 


5 : 
R= 3 R; Nu: (58) 
LL 


where N are the one dimensional shape functions (given in 
Appendix A). 
If we substitute Eq. 58 into 57 and make the coordinate 


transformation, after integration we get 


2 
2 


2 


of Ry 


zs rs OO me 2 4 “ 
Q,* Vat [a9 (42 2) (31R,+64R *46R,R,-8R,R 14R,R3). (59) 


I'2 13 


Now by equating Eq. 59 to Eq. 56 at any time tT we may deter- 
mine whether the front of the flow has already passed the 
mid-node of the wetted side of the element itl. 

A similar argument can be carried out for the end node 
of the element iti. In that case we also get an equation 
Similar to Eq. 59 with different numerical coefficients. 


For the case of the outside xlow:it has also been 


assumed to have a constant fluid flux and a fictitious 


entrance circular cross-sectional area (whose radius is the 
largest R plus unity) has been assumed. 

The numerical coefficients of R's in Eq. 59 are used 
in subroutine FLOW and it has been assumed that the temper- 
ature of a fluid "particle" does not change during passage 
through the active section. This is believed to be an 
acceptable approximation for representative values of flow 
velocity and active length. With this assumption, for a 
given entry time-temperature relation (inside and outside 
flow) this subroutine evaluates the fluid nodal temperatures 


at any time. 


6. SUBR@UTINE }kyRMV 

At every step of time integration this subroutine is 
called by the main program to form the vector v of the 
right-hand side of the discretized finite elemer. {.4. 13) 
for the given thermal boundary conditions The program 


itself is self explanatory. 


7.  SUBR@UTINE TEMPER 

The nodal temperatures are calculated in this subroutine 
with the trapezoidal time integrations. Irons' correction 
is applied here after every 10 steps of time integration. 
The temperature vectors selected for the stress analysis are 
stored. The transient temperatures will be printed out at 


the desired interval of time. 
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8. SUBROUTINE STIFF 

For any stress or thermal stress problem subroutine 
STIFF is called once by the main program to evaluate the 
stiffness matrix of the system. Once the stiffness matrix 
of the system is calculated then the desired structural 
boundary conditions are applied and, at the end, the Cholesky 
decomposition is performed. The functional flow chart of 


subroutine STIFF is given in Fig. 24. 


9. FORMF 

In subroutine F@RMF the thermal load vectors are calcu- 
lated for as many as (IVEC) given temperature vectors. 
These thermal load vectors are the columns of a rectangular 
matrix F. The provision is made that no matter how many 
temperature vectors are given (always IVEC < 20) the IVEC 
number of columns of the F are filled with the thermal load 
vectors and the last column of F will contain the load 
vector corresponding to the unit end displacement for zero 
axial force when the plane-end boundary condition is applied. 


The flow chart of subroutine F@RMF is given in Fig. 25. 


10. SUBR@UTINE CENTF 

If the system is rotating about the axis of revolution, 
this subroutine is called once by the main program to eval- 
uate the centrifugal load vector. This load vector is 
always placed in the first column after the thermal load 


vectors (IVEC+1 position) in F. 
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Figure 24, Functional Flow Chart o1 STIFF. 
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Figure 25. Functional Flow Chart of FORMF. 
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11. SUBR@UTINE PRESS 

If there is any pressure acting on the system, the 
pressure load vector is calculated in this subroutine and 
this vector is placed in the column next to centrifugal load 
vector, if any, otherwise it will fill the position desig- 


nated for centrifugal load vector in EF. 


12. SUBRGUTINE DISPL 

The displacement vectors are found in this subroutine. 
Since the stiffness matrix is already decomposed in banded 
form, then the displacement vectors one after another are 
obtained by back and forward substitution. The principle 
of superposition is applied here and finally the axial 
force, if any, is corrected for the plane-end boundary 


condition. 


13. SUBROUTINE STRESS 
For every problem this subroutine is called by the main 
program once to evaluate the stresses at the points 


(ntl, g= +4 


y)o£ each element. The transient stresses are 
calculated and printed for each displacement vector. For 
each problem the maximum mean stress and the maximum octa- 
hedral shearing stress, together with the corresponding 


times and locations, are printed. 


14, LIMITATIONS 
In the process of the development of the program dis- 
cussed so far, it has been intended that all of the calcula- 


tions and storage of data occur in the computer without 
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using any external devices such as disks or magnetic tapes 
so as to be able to solve any sizabis problem in relatively 
‘ short time. 
For this reason the following limitations (Table VI) are 
set forth which give an overall size of 450K bytes to the 


program. 


TABLE VI 
MAXIMUM VALUES FOR PROGRAM PARAMETERS 


NBAND Half band-width of the system stiffness matrix 66 


NMAT Total number of different materials 5 
NET Total number of elements 40 
NNT Total number of nodes 149 
’ NCF6T Total number of nodes in contact with 37 
outside fluid ; 
1 NCFIT Total number of nodes in contact with 37 
inside fluid 
NNRE Total number of nodes of right-hand end 37 
NNLT Total number of nodes constrained against any 37 
longitudinal motion 
NNRT Total number of nodes constrained against any 37 
radial motion 
IVEC Total number of temperature vectors stored for 20 
thermal stresses 
NRAMPG Total number of ramps for outside flow 15 
‘ NRAMPI Total number of ramps for inside flow 15 
NPNT Total number of pressure nodes af 
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